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We apply the recently developed critical minimum energy subspace scheme for the investigation 
of the random-field Ising model. We point out that this method is well suited for the study of 
this model. The density of states is obtained via the Wang-Landau and broad histogram methods 
in a unified implementation by employing the N-fold version of the Wang-Landau scheme. The 
random- fields are obtained from a bimodal distribution (hi = ±2), and the scaling of the specific 
heat maxima is studied on cubic lattices with sizes ranging from L — 4 to L = 32. Observing the 
finite-size scaling behavior of the maxima of the specific heats we examine the question of saturation 
of the specific heat. The lack of self-averaging of this quantity is fully illustrated and it is shown 
that this property may be related to the question mentioned above. 
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I. INTRODUCTION 

The random-field Ising model (RFIMMMis one of the 
best studied glassy magnetic models 0, 00 , mainly be- 
cause of its interest as a simple frustrated system. In fact, 
it has been a matter of conspicuous controversy over the 
last 15 years, mainly concerning the nature of its phase 
transitions. The RFIM Hamiltonian is given by: 

H = -J SiSj-J^hiSi (!) 

<i,j> i 

where Si = ±1, J > is the nearest- neighbors ferro- 
magnetic interaction and the random-fields (RF's) are 
obtained from a discrete distribution hi = ±A, where 
A = 2 is the disorder strength, also called randomness of 
the system. 

The notion of dimensional reduction [j| indicated that 
the critical behavior of the RFIM in d dimensions, at 
sufficiently low randomness, should be identical to that 
of the well-known normal Ising model in d— 2 dimensions. 
On the other hand, the droplet theory of domain wall 
energies in the ferromagnetic state |6( suggested that a 
phase transition should exist in three-dimensions (3D), 
for finite-temperature and randomness. The puzzle has 
been cleared out by Imbrie , Schwartz Q and Bricmont 
and Kupiainen 9] . Their arguments strongly support the 
view that a phase transition in 3D exists for sufficiently 
small randomness (A c w 2.3) 0. 

From the experimental point of view, a true realiza- 
tion of the RFIM is hardly conceived. However, it has 
been shown that dilute antiferromagnets in uniform ex- 
ternal field (DAFF) represent physical realizations of the 
RFIM 10] and a number of experiments investigated the 
phase transitions of such 3D systems ^l| . These experi- 
ments have proven to be very difficult and their interpre- 
tation doubtful due to the slow, glassy dynamics of the 
system. 

Although there exist several open questions about the 
phase transition in the RFIM, it is now generally ac- 



cepted that a new fixed point controls the behavior of 
RF ferromagnets 0, ^| . The significance of this for the 
RFIM (in d > 2) is that this new zero-temperature ran- 
dom fixed point controls the whole critical line (T C (A)) 
and that the RF's are always relevant. For disordered 
systems with weak randomness which couples to the lo- 
cal energy (such as random-site impurity or random-bond 
models) the crossover to a new random fixed point, de- 
pends on the Harris criterion [l3lll4| . According to this, 
the disorder is relevant if the correlation length expo- 
nent of the pure model (y = v pU re) satisfies the condition 
dv < 2 and this condition may be stated, with the help 
of the hyperscaling relation (a = 2 — dv), as a > 0. 
Since the specific heat exponent of the 3D Ising model 
is positive [l5j, weak disorder should be expected to be 
relevant. In the case of the RFIM the type of disorder is 
much more severe, since the randomness couples to the 
local order parameter and the crossover renormalization 
group eigenvalue is always positive fl3. The inequal- 
ity v > 2/d, derived by Chayes et al. for the cor- 
relation length exponent of a generic disordered system 
[y = Vrandom) would imply, using again hyperscaling, a 
negative specific heat exponent (a < 0). However, it is 
believed that hyperscaling is violated in the RFIM and 
the specific heat exponent a is related to v by a modi- 
fied hyperscaling law 2 — a = [d — &)v. The exponent 
9 characterizes the scaling of the stiffness of the ordered 
phase at the critical point |T^. Thus, the specific heat 
exponent of the RFIM is not restricted, by the above 
theoretical considerations, to be negative [12j . 

The inconsistency of various estimations in the litera- 
ture concerning the critical exponent a is the origin of a 
long lasting lively controversy, leaving open, so far, even 
the question of divergence or saturation of the specific 
heat. The specific heat of the RFIM can be experimen- 
tally measured and is of considerable theoretical interest. 
Several Monte Carlo methods at finite temperatures but 
also methods using ground state configurations have been 
used to estimate the critical exponent a. Some of the 
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ground state studies came u p w ith strongly negative val- 
ues, ranging from a = -1.5 |17| to a = -0.5 [la ITl l20| . 
whereas Middleton and Fisher estimated in marked 
disagreement a = —0.01 ± 0.09. Experiments on DAFF 
provided evidence of a second order phase transition and 
a logarithmic singularity for the specific heat [2jJ. Re- 
cently, Barber and Belanger in their Monte Carlo 
study of a DAFF model reported also that their specific 
heat curve closely mimics a logarithmic peak. Moreover, 
it has been pointed out that a strongly negative value 
of a causes serious difficulties when it comes to finding 
a consistent set of scaling relations to describe the crit- 
ical behavior of the RFIM. These scaling relations are 
consistent if one uses a sa [l^, which is also close to 
the experimental value |ll| . Clearly more work is needed 
to understand the specific heats behavior of the model. 
This important issue may be intimately linked to the 
main physical finding of this paper: the violation of self- 
averaging of the specific heat, illustrated below in Sec. lIIII 

The RFIM has been studied numerically using tradi- 
tional 0, 0, |U but also more sophisticated Monte 
Carlo techniques However, the nature of the model 
demands enormous computer resources. The equilibra- 
tion of the system at low temperatures is exponentially 
slow for large systems. Furthermore, in order to get a 
good estimate of the mean properties of the system, it 
is necessary to repeat the simulations for a large number 
of realizations of RF's. In the present work we have ap- 
plied the new and popular Wang-Landau (WL) [25| and 
broad histogram (BH) |26( methods to estimate the den- 
sity of states (DOS), G(E), of the model. These methods 
have been employed in a unified implementation u sing 
the Schulz et al. N-fold version of the WL scheme [27|, l28| 
and the energy space was restricted using the recent crit- 
ical minimum energy subspace (CMES) technique [l5| . 

The rest of the paper is organized as follows. In SeclLTI 
we provide an outline of the numerical methods which are 
involved in our calculations, including a brief description 
of the WL and BH methods. The recently developed 
CMES restriction is properly adapted and illustrated for 
the RFIM and useful technical details are provided. In 
Sec. IIIII we discuss the main conclusion of our work: the 
violation of self-averaging of the specific heat of the RFIM 
by studying the relevant probability distributions. The 
scaling behavior of the pseudocritical temperatures and 
their sample-to-sample fluctuations are also presented. 
Our conclusions are summarized in Sec. IIVI 



II. A NUMERICAL APPROACH 

We proceed to describe and appropriately adapt to the 
RFIM a recently developed Monte Carlo approach, based 
on the WL algorithm for estimating the DOS and us- 
ing the idea of dominant energy subspaces (CMES tech- 
nique). Consider a particular RF realization. Then, the 
specific heat and its peak are easily obtained with the 
help of the usual statistical sums. The CMES scheme 0] 



uses only a small but dominant part (E_, E + ) of the en- 
ergy space (E m i n , E max ) to determine the specific heat 
peaks. Let E denote the value of energy producing the 
maximum term in the partition function at the pseudo- 
critical temperature (corresponding to the specific heat 
peak) and S(E) = InG(E) the microcanonical entropy. 
Then, Eq. © defines the CMES approximation: 



C L (E-,E+) = N~ 1 T~ 2 iz-^J^exp [$(#)]- 

Z-^Eexp [$(£)] \ (2a) 



= [S(E)-f3E]- S(E) ~{3E , Z = ^ exp [$(£)] 

( 2b ) 

where N = L 3 and (E- , E+ ) is the minimum dominant 
subspace satisfying the following accuracy criterion: 



C L (E_ 1 E + ) 



C L (E, 



ni > n i ^max ) 



- 1 



< r 



(3) 



with r = 10 -6 . This accuracy is extremely demanding 
compared to the statistical errors produced by the DOS 
method and to the large sample-to-sample fluctuations 
of the RFIM. An algorithmic approach for specifying the 
CMES is described in Ref. 0. 

Using an ensemble of M(m = 1,...,M) macroscopic 
samples of linear size L corresponding to different RF 
realizations we have applied the described scheme in a 
broad energy space that covers the overlap of the domi- 
nant energy subspaces for all RF's of the ensemble. For 
a RF realization, say m, let us denote by (EL. m , E +>m ) 
the location of the dominant energy subspace defined 
by the above restriction, and by AE m — E +t . m — E-. m 
its extension. Our simulations were carried out in the 
broad energy subspace, which covers at least the union 
of the individual subspaces, i.e. (EL[m], -E'+rMl) = 
U m (J5_, m , E+, m ), of total extension (AE) [M] = S+JM] - 
E_ [m]- Note that, the total extension may be for large 
lattices several times larger than the individual exten- 
sions. This practice has the advantage that the approxi- 
mation of the specific heat for a particular RF is accurate 
in a wide temperature range, including its pseudocrit- 
ical temperature. Thus, the present implementation is 
not the most efficient for the purposes of locating only 
the specific heat peaks. However, this usage provides a 
more reliable alternative for comparing the statistics of 
the specific heat peaks to the averaged specific heat curve 
used in the literature (see for example Ref. 0] and also 
Sec. IHIfl and for discussing the pathology of this quite 
common choice. Despite the strong fluctuations of the 
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energy value corresponding to the maximum term of the 
partition function Z, the union of the CMES for large 
samples of RF's is a relatively small subspace, compared 
to that of the normal Ising model. Consider the case 
L = 16. Then, the energy levels used in Ref. for 
the normal Ising model, counting from the ground state 
(ie = 1), are the levels (ie = 1220 - 2410). Using an 
ensemble of 1000 RF's the union of the CMES was found 
to be the range of levels (ie = 1 — 950), while our simu- 
lations were performed in a wider range (ie = 1 — 1200). 
This is of the same order with that of the normal Ising 
model and at least five times smaller than that of the 
total energy space. 
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FIG. 1: The specific heat for two characteristic examples of 
RF's, RF a (dotted line) and RF b (solid line). The specific 
heat curves were obtained using the average DOS over the 
50 runs for each RF. The inset shows the fluctuation of the 
extension of the individual dominant subspaces over 50 runs. 

To conclude the above technical remarks, let us illus- 
trate that the efficiency of our method may be highly 
increased by studying each RF realization in its own dom- 
inant subspace. The following observations will be useful 
in subsequent studies of the RFIM or analogous mod- 
els. Fig. shows results of an extensive WL simulation 
of two particular RF's, labelled as RF m=a = RF a and 
RF m= b = RF\j in the figure, for a lattice of linear size 
F = 24. The simulation was repeated, for both RF's, 
50 times in the energy subspace (ie — 200 — 2000). The 
curves of the specific heat shown were obtained from the 
average DOS over the 50 runs. Note however that, the 
union space used in our simulations presented in Sec. IIIII 
(Fig. |3J) was estimated over large ensembles of RF's. For 
instance, the union space for F = 24 and its exten- 
sion was found to be of the order of 2500 energy lev- 
els (ie = 200 - 2700), using an ensemble of 200 RF's. 
Since there are RF's with one sharp peak and RF's with 
two or more pronounced peaks we have chosen to show in 
Fig.Q]two characteristic examples of RF's, RF a and RFb- 
The inset of Fig. ^ shows the fluctuation of the extension 
of the individual dominant subspaces over the 50 runs. 



There are some points that one should observe from this 
figure. Firstly, in both cases the extension of the CMES 
for a particular RF is much smaller than the broad en- 
ergy space used in the simulations. Specifically, for the 
RF a the dominant energy subspace is approximately of 
the order of 850 energy levels (ie = 800 — 1650), which 
is almost three times smaller than the 2500 levels of the 
total union space. Secondly, the fluctuations of the ex- 
tension AEb of the RFb are more pronounced and this is 
related to the existence of a secondary peak in the left 
of the main peak. This secondary peak causes a stronger 
fluctuation in the estimation of the end points of the cor- 
responding dominant subspaces. In any case, we could 
improve the efficiency of our scheme by a factor of at least 
two (for the case F — 24), by carefully individualizing the 
used energy space for simulating a particular RF. In fact, 
the fast early stages of the WL process may be used as 
a prognostic method to approximately locate the CMES 
of a particular RF, and this strategy may be an indis- 
pensable ingredient in analogous future studies. Finally, 
an entropic sampling study of the magnetic properties of 
the RFIM using the CMES restrictive entropic scheme 
based on the high-levels of the WL algorithm |2{| would 
be greatly facilitated by such a strategy. 

To determine the density of states, we have used the N- 
fold version of the WL method as presented by Schulz et 
al. [23, using 20 iterations for the reduction (fj+i = 
fi — e ) °f the WL modification factor. Our im- 
plementation is analogous to that presented in Ref. [28| . 
where the first 13 iterations follow the simple WL scheme 
and the rest iterations (j = 14—20) use the N-fold version 
of Schulz et al. [27| ■ We have used a flatness criterion of 
5% for the energy histogram [IH^l]. For RF's hi = ±2 
the classes for the N-fold process are specified by the en- 
ergy changes AE n = -16 + 4 • (n — l),n= 1, 2, 9. Us- 
ing the part of the simulation corresponding to the N-fold 
iterations (j = 14 — 20), we have accumulated data corre- 
sponding to non-zero energy changes in order to apply the 
well-known BH equation j26|: G(E)(N(E, E+AE n )) E = 
G(E + AE n )(N(E + AE n , E)) E+AEn . N(E, E + AE n ) is 
the number of possible spin flip moves from a microstate 
of energy E to a macrostate with energy E + AE n , which 
are known during the N-fold process. In this way we have 
produced 4 BH (n — 1,2,3,4) approximations for the 
DOS and the specific heat for each RF of the ensemble. 



III. LACK OF SELF-AVERAGING OF THE 
SPECIFIC HEAT 

For a disordered system we have to perform two dis- 
tinct kinds of averaging. For each sample, the usual ther- 
mal average has to be carried out and then we have to 
average over the random parameters. Let C m (T) de- 
note the specific heat of a particular realization m in the 
ensemble of M realizations of RF's. The pseudocritical 
temperature T£ m will, of course, depend on the realiza- 
tion of the RF. The location of the corresponding peak 
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is denoted by (C^,T£ m ) and the respective probability 
distributions by P L (C^) and P L (T* L J. 
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FIG. 2: Averaged curves for various samples of RF's. 

Rieger and Young 0, 0] have studied the following 
sample summation for the specific heat curves: 



1 M 

[C]av = -JJ ^ C m (T) 
m—1 



(4) 



and the finite-size scaling behavior of the peak of this 
averaged curve has been studied by assuming that the 
maximum [C]* v — maxT[C] av and the corresponding 
pseudocritical temperature T£ obey the scaling laws: 



T r + bL 



-1/v 



(5a) 



(5b) 



curve develops, as we increase the number of RF's. 
Fig. |3 presents an example of the probability distribu- 
tion PL(T£ m ). Since the peaks are found in different 
locations, the averaging in (0J wipes the particular peaks 
out. This explains why the averaged curve does not rep- 
resent the behavior of the most probable, say x, realiza- 
tion of the RF's: [C] av (T) ^ C X (T). It also suggests 
the absence of self-averaging for the specific heat of the 
present model, at least for the randomness studied here. 
Fig. Q] illustrates the finite-size behavior of the distribu- 
tion Pi(C* ra ). Although for L = 4 the distribution is 
sharp, as L increases the distribution broadens so that 
there is a significant number of RF's having their peaks 
higher, or lower, than the expected sample mean, defined 
in Eq. (|6al) . The above observations provide very strong 
evidence that the real behavior of the RFIM is not ap- 
propriately described by a possible misleading saturation 
of [C]*„. The source of this problem is the severe fluctu- 
ation of the pseudocritical temperatures and the lack of 
self- averaging may be an important statement. 

Broad distributions, with lack of self- averaging have 
been studied also in other physical problems, such as in 
the well-known case of the scaling theory of Anderson 
localization. There, it has been shown that for a dis- 
order electronic sample the conductance distribution at 
the point of the metal-insulator transition (the mobil- 
ity edge) is so broad, that the conductance is not a self- 
averaging quantity |3(j . Noteworthy that, the lack of self- 
averaging appears to be a common property of disordered 
systems at criticality and that besides the above men- 
tioned paradigm, one can find several ex amp les of mag- 
netic systems where this feature is present [3lll32ll33 . l34j . 
Actually, when dealing with physical quantities that are 
characterized by broad distributions, one must be mind- 
ful when attempting to define a transition in terms of 
related averaged quantities. In this sense, it seems that 
for the present model the common use of [C']* v may be a 
"meaningless" choice for a proper description. In order 



Here we shall also examine the scaling of the sample av- 
erages of the specific heat maxima and the pseudocritical 
temperatures, defined by: 

[C* m U = Yt £ C *m = P + cL^ ( 6a) 



L,m]av = TT ^ ^L,m — T c + bL X I V 



(6b) 



The possibility of different exponents may be ulti- 
mately related to the functional form of the distributions 
Pl(CJ,) and Pl(T£ m ), whose behavior is decisive for the 
comprehension of the critical behavior of the RFIM. 

For a small number of RF's the averaged curve [C] av 
has several local maxima reflecting a very strong sample- 
to-sample fluctuation of the individual pseudocritical 
temperature. Fig. [21 shows how the smoothness of this 
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FIG. 3: Fluctuation of pseudocritical temperature, 
to discuss the significance of the above broad probability 
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distributions and to present a more convincing finite-size 
scaling argument for the violation of self-averaging in the 
thermodynamic limit, we have included as an inset in 
Fig. |I] the ratio R c = V C /[C^^ V , where V c is the samplc- 
to-sample variance of the average (|Ha)l . Both the WL and 
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FIG. 4: Broadening of probability distributions. The inset 
presents finite-size evidence for the violation of self-averaging. 
R c is denned and discussed in the text. 



specific heat for most RF's. Nevertheless, this observed 
underestimation was not of the order of the decline in 
Fig. [SI so it is possible that the model crossovers to the 
conjectured saturation at these lattice sizes. 

From our attempts to acquire a better comprehension 
of the reasons for the above mentioned underestimation 
we also observed that this aspect is quite strong for RF's 
with a sharp specific heat peak. Note that such RF's 
are quite common and have been recently discussed also 
by Wu and Machta 37]. This underestimation may be 
observed also within the j = 20 WL level by using multi- 
ple measurements and also a separation S — 16 between 
successive recordings of the accepted microstates of the 
WL process. The separation refinement is generally be- 
lieved to improve the accuracy of the WL method HlH. 
Fig. [5] illustrates its effect in a repeated application us- 
ing the RF a , that appears also in Fig. 2] From Fig. 
we observe that the effect of separation is to increase 
the mean value of the maximum of the specific heat by 
an amount which is of the same order with the stan- 
dard deviation of the statistical errors that one obtains 
by using multiple measurements (100 independent WL 
runs) without separation. The standard deviation of the 



BH estimates are shown with their errors. The variance 
V c was reduced, by eliminating the statistical (DOS) er- 
rors, assumed to be of the order of the difference between 
the two methods. The above defined normalized square 
width is a measure characterizing the self-averaging prop- 
erty of a system [33, E3 HH • This ratio appears to tend 
to a constant (R c — ► 0.3), as can be seen from the inset of 
Fig. H Thus, according to the literature [H the 
system is not self-averaging and the corresponding dis- 
tribution does not become sharp in the thermodynamic 
limit. 

The WL and BH estimates for [C^Jou appear in Fig.0 
For L = 4 — 20 we have averaged over an ensemble of 1000 
RF's and for L > 20 over 200 RF's. In this figure we show 
the WL estimates and the mean (BH) of the 4 BH esti- 
mates. The same figure presents the size dependence of 
[C]* v . Although for the range L = 4 - 20 the behavior of 
the estimates is convincing for their accuracy, an increase 
of statistical errors is observed for larger sizes, depicted 
in the growing differences between the WL and BH esti- 
mates. The estimates for L = 28 and L — 32 appear to 
decline from the L = 4 — 20 behavior and the growing 
errors after L — 24 make difficult a definite judgment for 
the asymptotic behavior. Refinements of the WL scheme 
will be favorable for these larger lattice sizes. This could 
be attempted by using multiple measurements for each 
RF, by increasing the final WL j-iteration and/or by in- 
troducing other refinements of the WL algorithm, such 
as a separation S between successive recordings |36] . Our 
first attempt to increase the WL j-iteration to j — 24 for 
small samples of RF's indicated that the level of j = 20 
leads to an underestimation of the sharp peaks of the 




FIG. 5: Size behavior of the averages [CJ,] m and The 
vertical bars illustrate the order of the sample-to-sample fluc- 
tuations and should not be confused with the small errors of 
the WL scheme. 

new sample of multiple measurements (25 independent 
WL runs using separation) is also of the same order, as 
shown in Fig. El We note that, the WL sampling in these 
multiple measurements was carried out in a energy sub- 
space which is slightly wider than the CMES of the RF a 
(ie = 700 — 1800) and not in the wider energy range 
used for the simulation appearing in Fig. ^ Comparing 
these two figures (Fig. ^ an d Fig. EJ) one can detect the 
effects of different restrictions on the energy space. The 
observed sadden decrease of the right tail of the specific 
heat (a similar comment applies also for the left tail) 
in Fig. ^ is an effect induced by the restriction imposed 
on the energy space and appears in the neighborhood of 
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T w 3. The further restriction imposed in the new sam- 
plings (appearing in Fig. |SJ is now reflected in the shift 
of the observed sadden decrease in the neighborhood of 
T Kt 2.5. Before attempting to simulate larger samples 
of RF's, other refinements should be also tested, in or- 
der to obtain a more accurate and efficient scheme. In 
any case, our study shows that there is a large number of 
RF's with sharp peaks strongly fluctuating in their pseu- 
docritical temperatures and this generic property makes 
the proposed CMES scheme the most appropriate alter- 
native, despite the accuracy and slowing down problems 
observed at the larger sizes. Sharp peaks are usually 
missed by importance sampling, due to an inadequate 
temperature scanning often used. 
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FIG. 6: Illustration of the separation effect in the specific heat 
of the RF a . The specific heat curves shown were calculated 
from the average DOS over the runs, while the horizontal lines 
represent the mean values of the independent runs for the 
specific heat peaks. The dotted line shows the case S = 16, 
while the solid line the case 5 = 0. The error bars illustrate 
the standard deviation (of the independent peaks) over the 25 
(dotted) and 100 (solid) WL runs, corresponding to S = 16 
and S — 0. 

Fig. [7| illustrates that the two pseudocritical tempera- 
tures, [T£ m ]av and T£, tend to the same limit (T c = T c ). 
The behavior of [T£ m ] av is smoother than the behavior 
of T£, which is more sensitive to the sample size. Us- 
ing our data for [T£ m ] a « we found a reasonably good fit 

with f c = 2.03(18) and v = 1.31(18). This value is very 
close to the value v — 1.37(9) found in Ref. [Ig and lies 
between the values 1.0(1) of Ref. |2(j and the estimates 
1.6(3) and 1.4(2) of Ref. [l7|. The inset in Fig. □ illus- 
trates the scaling of the sample-to-sample variance of the 



average [T£ m ] av of Eq. JBdJ). Assuming that the square 
of these sample-to-sample fluctuations scales with the lin- 
ear size L according to <5 2 ([T£ m ] a v) ~ L~ 2 / u , we obtain, 
from the fit shown in the inset, the value v = 1.18(15). 
This estimate is slightly smaller than the value found 
above, which is in good agreement with the best esti- 



mate in literature 0] . The rather slow approach of the 
fluctuations to zero is also an interesting finding. Ac- 
cording to Aharony and Harris [33| and Wiseman and 
Domany [34| , the fact that the square width of the dis- 
tribution of the sample dependent pseudocritical temper- 
atures scales with L~ 2 / v and not with L~ d , when com- 
bined with finite-size scaling |34| , is an indication of lack 
of self-averaging of the random system. Therefore, our 
main conclusion is reinforced and is also in conformity 
with the results of Parisi and Sourlas [3^| . The numeri- 
cal study of these authors showed that the strong fluctu- 
ations of the 3D RFIM produce a maximal violation of 
self- averaging for the correlation length. It appears that 
the disorder present in the RFIM brings about drastic ef- 
fects and its strong non self-averaging behavior includes 
also the specific heat, as suggested in this paper. 

It is quite possible that the above relevant aspect was 
overlooked in previous finite-temperature studies but also 
in ground state calculations. We think that, at least 
partly, this practice is behind the existing controversial 
situation in the literature concerning the behavior of the 
specific heat. For instance, the saturation of \C]* av is 
evident from the very small sizes and its behavior does 
not admit a finite-size scaling, but rather appears as a 
random fluctuation around the value 0.84, as shown by 
the dashed line in Fig. [SJ The negative value for the 
exponent a found, from the study of [C]* v in previous 
finite-temperature studies 0, 0] seems to us question- 
able. In addition to all the reasons mentioned above, 
the very early and clear saturation observed here and 
the possibility of a crossover behavior of the model at 
larger lattice sizes are strong indications that make us 
question the meaning of such a scaling prediction. It ap- 
pears that the behavior of [C^J a „ incorporates more of 
the physical content of the model, although its asymp- 
totic behavior seems unsettled at the lattice sizes stud- 
ied here. As pointed out earlier, the systematic errors 




FIG. 7: Size dependence of pseudocritical temperatures. The 
vertical bars as in Fig. Q>] The inset illustrates the scal- 
ing of the sample-to-sample variance of the average [T£ m ] u 
(Eq. (HP). 
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of the WL scheme for large lattice sizes are due to the 
practical (j — 20) but not fully converged usage of this 
algorithm in our simulations. This option was dictated 
by the need of studying large samples of RF's. We as- 
sume that the decline of the estimates observed here for 
L > 24 is stronger, as pointed out earlier, from the un- 
derestimation observed by studying smaller samples of 
RF's and using longer (j = 24) runs. Then, it is quite 
obvious from Fig. |3] that the true asymptotic behavior 
cannot be observed at these lattice sizes, although its 
saturation seems to be now plausible. In order to ob- 
tain a safe and sound estimation of the large L-behavior, 
larger systems of at least of the order of L = 60 should 
be considered. This is an extremely demanding computer 
project and will have to be postponed, until further tests 
make available a more accurate and optimum refinement 
of the presented scheme. Finally, let us point out that, 
our first attempts to observe the behavior of the suscep- 
tibility of the 3D RFIM via a recently proposed entropic 
scheme [2^|, suggested also an even stronger violation of 
self-averaging for the magnetic properties of the system. 



of the specific heat of the RFIM is still an open impor- 
tant topic, necessary for a better comprehension of the 
model. This problem was considered in a completely new 
basis in this paper. The property of self-averaging of 
the specific heat was addressed in a concise way and its 
violation was explicitly shown by studying the relevant 
probability distributions. This finding may lead to a bet- 
ter theoretical and numerical approach of the problem. 
The scaling behavior of the pseudocritical temperatures 
and their sample-to-sample fluctuations were also pre- 
sented, and found to support a strong violation of the 
self- averaging property of the system. The new ideas 
and numerical techniques utilized to tackle the RFIM 
use as an essential ingredient the critical minimum en- 
ergy subspace scheme. We hope that the combination 
of algorithms and techniques applied here will be useful 
in further numerical studies of this and other similarly 
challenging problems. 



IV. CONCLUSIONS 
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